package com.asa.numanaly;

import com.asa.HanShu;

/**
 * 插值和多项式逼近
 * p92
 * @author Administrator
 *
 */
public class C {

	
	/**
	 * Neville迭代差值法
	 * @param x
	 * @param y
	 * @param xpoint
	 * @return	结果是只需要取最后一排的
	 */
	public static Double[][] chazhi1(double[] x,double[] y,double xpoint) {
		
		Double[][] result = new Double[x.length][x.length];
		
		for (int i = 0; i < y.length; i++) {
			result[i][0] = y[i];
		}
		
		
		
		for (int i = 1; i < result.length; i++) {
			for (int j = 1; j < i; j++) {
//				
				result[i][j] = ((xpoint-x[i-j])*result[i][j-1]-(xpoint-x[i])*result[i-1][j-1])/(x[i]-x[i-j]);
				//System.out.print(result[i][j]);
			}
		}
		
		
		return result;
	}
	
	/**
	 * 将chazhi1转化为可以运行的函数
	 * @param c	chazhi1返回参数
	 * @return	函数
	 */
	public static HanShu chazhi12hanshu(final  Double[][] c){
		
		HanShu result = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x) {
				// TODO Auto-generated method stub
				double asa = 0;
				for (int i = 0; i < c[c.length-1].length-1; i++) {
//					System.out.println(c[c.length-1][i] );
					asa += c[c.length-1][i] * Math.pow(x, i);
				}
				
				return asa;
			}
		};
		
		return result;
		
		
		
	}

	
	
	
	
	
	
	
	
	
	
	
	/**
	 * 差商型Newton 插值公式	*
	 * p110
	 * @param x
	 * @param y
	 * @param xpoint
	 * @return
	 */
	public static Double[][] chazhi2(double[] x,double[] y,double xpoint) {
		
		
		Double[][] result = new Double[x.length][x.length];
		
		for (int i = 0; i < y.length; i++) {
			result[i][0] = y[i];
		}
		
		
		
		for (int i = 1; i < result.length; i++) {
			for (int j = 1; j < i; j++) {	
				result[i][j] = (result[i][j-1]-result[i-1][j-1])/(x[i]-x[i-j]);
			}
		}
		
		return result;

	}

	/**
	 * Hermite插值	*
	 * 
	 * @param x
	 * @param y
	 * @param ydao y这一点的导数,这一点的导数是让人迷惑的地方，这个数列的原函数我们是未知的，那么这一点的导数如何求呢？  例子中都是直接给出来相应点的导数值得
	 * 可能是这种方法固有的限制吧，好奇怪，之前我想过去用两个相邻点进行求两点斜率的方式来拟合求导（这个方法就是正确的方法，我们求一个函数在某
	 * 点的导数时，我们并不需要知道这个函数的样子，只需要知道一系列的点），但书中的东西和我这个观点又有出入
	 * 反正还有别的拟合法，比如三次样条插值没这限制，如果硬要用这种方法，而又没有给出相应点的导数（这在实际问题上是很可能出现的），
	 * 那就按我想的来吧
	 * 
	 * @return Q00，Q11。。。Q(2n+1,2n+1)
	 * 
	 * H(x) = Q00 +Q11(x-x0)+Q22(x-x0)^2+Q33(x-x0)^2(x-x1)
	 * 	+Q44(x-x0)^2(x-x1)^2+...
	 * 	+Q[2n+1,2n+1](x-x0)^2(x-x1)^2(x-x2)^2...(x-x[n-1])^(x-xn)
	 * 
	 */
	public static Double[] chazhi3(double[] x,double[] y,double[] ydao) {
		
		Double[] z = new Double[x.length*2+1];
		Double[][] result = new Double[x.length*2+1][x.length*2+1];
		
//		for (int i = 0; i < y.length; i++) {
//			result[i][0] = y[i];
//		}
		
		
		
		for (int i = 0; i < result.length; i++) {

			z[2*i] = x[i];
			z[2*i+1] = x[i];
			result[2*i][0] = y[i];
			result[2*i+1][0] = y[i];
			result[2*i+1][1] = ydao[i];

			if (i!=0) {
				result[2*i][1] = (result[2*i][0]-result[2*i-1][0])/(z[2*i]-z[2*i-1]);
			}
			
		}
		for (int i = 2; i < z.length; i++) {
			for (int j = 2; j < i; j++) {
				result[i][j] = (result[i][j-1]-result[i-1][j-1])/(z[i]-z[i-j]);

			}
		}
		Double[] Q = new Double[x.length*2+1];

		for (int i = 0; i < Q.length; i++) {
			Q[i] = result[i][i];
		}
		
		return Q;

	}
	
	
	
	/**
	 * 自然三次样条	*
	 * p129
	 * @param x
	 * @param y
	 * @param xpoint
	 * @return S(x)方程的四个参数abcd
	 * 这个的计算最终结果是一段又一段的分段函数拼凑起来的
	 * 	其中S(x) = S[j](x) = a[j]+b[j](x-x[j])+c[j]*(x-x[j])^2+d[j]*(x-x[j])^3	 其中x[j]<x<x[j+1]	j最大是n-1
	 */
	public static Double[][] chazhi4(double[] x,double[] y) {
		
		double[] h = new double[x.length-1];
		Double[] a = new Double[x.length-1];
		
		double[] l = new double[x.length];
		double[] u = new double[x.length];
		double[] z = new double[x.length];

		for (int i = 0; i < x.length-1; i++) {
			
			h[i] = x[i+1]-x[i];
			
		}
		
		for (int i = 0; i < x.length-1; i++) {
			a[i] = 3*(y[i+1]-y[i])/h[i] - 3*(y[i]-y[i-1])/3;
		}
		l[0] = 1;
		u[0] = 0;
		z[0] = 0;
		
		for (int i = 1; i < x.length-1; i++) {
			
			l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*u[i-1];
			u[i] = h[i]/l[i];
			z[i] = (a[i]-h[i-1]*z[i-1])/l[i];
			
		}
		l[l.length-1] = 1;
		z[z.length-1] = 0;
		Double[] c = new Double[x.length];

		c[c.length-1] = 0.0;
		
		Double[] b = new Double[x.length];
		Double[] d = new Double[x.length];

		for (int j = x.length-2; j >=0; j--) {
			c[j] = z[j]-u[j]*c[j+1];
			b[j] = (y[j+1]-y[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3;
			d[j] = (c[j+1]-c[j])/(3*h[j]);
		}
		Double[][] result = new Double[4][]; 
		result[0] = a;
		result[1] = b;
		result[2] = c;
		result[3] = d;
		return result;
		
	}
	
	
	
	
	public static HanShu chazhi42hanshu(final  Double[][] c){
		
		HanShu result = new HanShu() {
			
			@Override
			public double hanshu(double t, double x) {
				// TODO Auto-generated method stub
				return 0;
			}
			
			@Override
			public double hanshu(double x) {
				// TODO Auto-generated method stub
				
				
				
//				if (condition) {
//					
//				}
				
				
				
				double asa = 0;
				for (int i = 0; i < c[c.length-1].length-1; i++) {
//					System.out.println(c[c.length-1][i] );
					asa += c[c.length-1][i] * Math.pow(x, i);
				}
				
				return asa;
			}
		};
		
		return result;
	}

	
	
	
	
	
	
	
	
	
	
	/**
	 * 固支三次样条	*
	 * p131
	 * @param x
	 * @param y
	 * @param FPO	f(x0)的导数
	 * @param FPN	f(xn)的导数
	 * @return S(x)方程的四个参数abcd
	 * 	其中S(x) = S[j](x) = a[j]+b[j](x-x[j])+c[j]*(x-x[j])^2+d[j]*(x-x[j])^3		j最大是n-1
	 */
	public static Double[][] chazhi4_2(double[] x,double[] y,double FPO,double FPN) {
		
		double[] h = new double[x.length-1];
		Double[] a = new Double[x.length];
		
		double[] l = new double[x.length];
		double[] u = new double[x.length];
		double[] z = new double[x.length];

		for (int i = 0; i < x.length-1; i++) {
			
			h[i] = x[i+1]-x[i];
			
		}
		
		a[0] = 3*(y[1]-y[0])/h[0] - 3*FPO;
		a[a.length-1] = 3*FPN - 3*(y[y.length-1]-a[y.length-2])/h[h.length-1];

		
		for (int i = 1; i < x.length-1; i++) {
			a[i] = 3*(a[i+1]-a[i])/h[i] - 3*(a[i]-a[i-1])/h[i-1];
		}
		l[0] = 2*h[0];
		u[0] = 0.5;
		z[0] = a[0]/l[0];
		
		for (int i = 1; i < x.length-1; i++) {
			
			l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*u[i-1];
			u[i] = h[i]/l[i];
			z[i] = (a[i]-h[i-1]*z[i-1])/l[i];
			
		}
		l[l.length-1] = h[h.length-1]*(2-u[u.length-1]);
		z[z.length-1] = (a[a.length-1]-h[h.length-1]*z[z.length-1])/l[l.length-1];
		Double[] c = new Double[x.length];

		c[c.length-1] = z[z.length-1];
		
		Double[] b = new Double[x.length-1];
		Double[] d = new Double[x.length-1];

		for (int j = x.length-2; j >=0; j--) {
			c[j] = z[j]-u[j]*c[j+1];
			b[j] = (a[j+1]-a[j])/h[j] - h[j]*(c[j+1]+2*c[j])/3;
			d[j] = (c[j+1]-c[j])/(3*h[j]);
		}
		Double[][] result = new Double[4][]; 
		result[0] = a;
		result[1] = b;
		result[2] = c;
		result[3] = d;
		return result;
		
	}
	
	/**
	 * Bezier曲线	*
	 * p141
	 * @param x
	 * @param y
	 * @param FPO	f(x0)的导数
	 * @param FPN	f(xn)的导数
	 * @return (xi(t),yi(t))方程的四个参数abcd
	 * 	其中(xi(t),yi(t)) = (a0[i]+a1[i]*t+a2[i]*t*t+a3[i]*t*t*t,b0[i]+b1[i]*t+b2[i]*t*t+b3[i]*t*t*t)	j最大是n-1
	 */
	
	public static Double[][] chazhi5(double[] x,double[] y,double[] xp,double[] yp,double[] xm,double[] ym) {
		double[] a0 = new double[x.length-1];
		double[] a1 = new double[x.length-1];
		double[] a2 = new double[x.length-1];
		double[] a3 = new double[x.length-1];
		double[] b0 = new double[x.length-1];
		double[] b1 = new double[x.length-1];
		double[] b2 = new double[x.length-1];
		double[] b3 = new double[x.length-1];
		
		for (int i = 0; i < a0.length; i++) {
			
			a0[i] = x[i];
			b0[i] = y[i];
			a1[i] = 3*(xp[i]-x[i]);
			b1[i] = 3*(yp[i]-y[i]);

			a2[i] =  3*(x[i]+xm[i+1]-2*xp[i]);
			b2[i] =  3*(y[i]+ym[i+1]-2*yp[i]);

			a3[i] = x[i+1]-x[i]+3*xp[i]-3*xm[i+1];
			b3[i] = y[i+1]-y[i]+3*yp[i]-3*ym[i+1];

		}
		
		
		return null;//没写完,返回所有a和b

	}

	
	
	
	
	
	
	
	
	
	
	
	
	
	
	public static void main(String[] args) {
		test1();
		
		
		
		
		
		
		
	}








	private static void test1() {
		double[] x = new double[]{8.1,8.3,8.6,8.7};
		double[] y = new double[]{16.99410,17.56492,18.50515,18.82091};
		Double[][] chazhi1 = chazhi1(x, y,8.4);
		
		
		for (int i = 0; i < chazhi1.length; i++) {
			for (int j = 0; j < chazhi1[i].length; j++) {
				System.out.print(chazhi1[i][j] +"      ");
			}
			System.out.println();
		}
		
		
		
	}
	
	
	
	
	
	
}
